Integrated analysis of the transcriptome, sRNAome, and degradome reveals the network regulating fruit skin coloration in sponge gourd (Luffa cylindrica)

Sponge gourd fruit skin color is an important quality-related trait because it substantially influences consumer preferences. However, little is known about the miRNAs and genes regulating sponge gourd fruit skin coloration. This study involved an integrated analysis of the transcriptome, sRNAome, and degradome of sponge gourd fruit skins with green skin (GS) and white skin (WS). A total of 4,331 genes were differentially expressed between the GS and WS, with 2,442 down-regulated and 1,889 up-regulated genes in WS. The crucial genes involved in chlorophyll metabolism, chloroplast development, and chloroplast protection were identified (e.g., HEMA, CHLM, CRD1, POR, CAO, CLH, SGR, CAB, BEL1-like, KNAT, ARF, and peroxidase genes). Additionally, 167 differentially expressed miRNAs were identified, with 70 up-regulated and 97 down-regulated miRNAs in WS. Degradome sequencing identified 125 differentially expressed miRNAs and their 521 differentially expressed target genes. The miR156, miR159, miR166, miR167, miR172, and miR393 targeted the genes involved in chlorophyll metabolism, chloroplast development, and chloroplast protection. Moreover, a flavonoid biosynthesis regulatory network was established involving miR159, miR166, miR169, miR319, miR390, miR396, and their targets CHS, 4CL, bHLH, and MYB. The qRT-PCR data for the differentially expressed genes were generally consistent with the transcriptome results. Subcellular localization analysis of selected proteins revealed their locations in different cellular compartments, including nucleus, cytoplasm and endoplasmic reticulum. The study findings revealed the important miRNAs, their target genes, and the regulatory network controlling fruit skin coloration in sponge gourd.

www.nature.com/scientificreports/ it considerably influences consumer preferences. To date, many genes related to fruit skin coloration have been reported for horticultural plants. For example, genes controlling the fruit skin colors of cucumber fruits, including white immature fruits (w and w 0 ) 10,11 , orange mature fruits (B) 12 , and fruits with a yellowish green peel (ygp) 13 , have been identified and mapped. Additionally, several genes controlling chloroplast development and chlorophyll biosynthesis have been detected in tomato fruits. The protein encoded by SlMYB72 directly targets genes involved in the metabolism of chlorophylls, carotenoids, and flavonoids, thereby determining tomato fruit color and ripening 14 . The GREEN STRIPE (GS) gene, which is a methylated isoform of TOMATO AGAMOUS-LIKE 1 (TAGL1), reportedly regulates diverse chloroplast developmental processes and carotenoid accumulation in tomato fruit 15 . Moreover, BEL1-LIKE HOMEODOMAIN 11 (SlBEL11) 16 , KNOTTED1-LIKE HOMEOBOX (KNOX) genes (TKN2 and TKN4) 17 , ARABIDOPSIS PSEUDO RESPONSE REGULATOR 2-LIKE (SlAPRR2-LIKE) 18 , and the auxin response factor (ARF) gene SlARF4 19,20 also affect fruit chloroplast development or chlorophyll biosynthesis in tomato. However, the critical genes involved in sponge gourd fruit skin coloration remain unknown. MicroRNAs (miRNAs) are a class of noncoding short RNAs (19-24 nt) that participate in post-transcriptional regulation by inhibiting translation or cleaving targeted messenger RNAs (mRNAs) 21 . Previous studies revealed a role for several miRNAs in the tissue coloration of horticultural plants [22][23][24][25][26] . Specifically, miR156 and its target SQUAMOSA PROMOTER BINDING PROTEIN-LIKE (SPL) genes contribute to the regulation of light-induced red peel coloration and anthocyanin accumulation in pear 22 . The MIR156a-SPL12 module coordinates the accumulation of chlorophylls and anthocyanins during blueberry fruit ripening 23 . An earlier study proved that miR828 regulates phenylpropanoid accumulation by modulating the expression of R2R3-MYB transcription factor-encoding genes and is associated with anthocyanin production in potato 24 . Another study demonstrated that miR858 negatively regulates anthocyanin biosynthesis by repressing the expression of AaMYBC1 in red kiwifruit 25 . Moreover, the overexpression of bol-miR171b in broccoli results in dark green leaves with increased chlorophyll contents 26 .
Combined multi-omics analyses enabled us to identify several crucial genes and their regulatory network associated with pigment formation in horticultural plants 27,28 . In cucumber, the chlorophyll biosynthetic genes associated with decreased chlorophyll or chloroplast levels in white fruit skin were identified and a predicted anthocyanin biosynthesis regulatory network was established on the basis of an integrated analysis of the metabolome and transcriptome 27 . De novo transcriptome and metabolome analyses of green and purple turnips detected anthocyanins and key genes mediating the difference in root skin pigmentation 28 . Transcriptome, sRNAome, and degradome analyses have been combined to investigate the crucial miRNAs and their target genes as well as the associated network in plants [29][30][31] . An integrated transcriptome, sRNAome, and degradome examination revealed that miRNAs can mediate tea plant immunity by regulating differentially expressed genes (DEGs) at the post-transcriptional level and miR530b-ERF96 (encoding ethylene response factor 96) and miRn211-TLP (encoding thaumatin-like protein) are important for the responses to gray blight 30 . The high-throughput sequencing of the sRNAome, degradome, and transcriptome identified a novel regulatory pathway involving ethylene-miR164-NAC that modulates kiwifruit ripening 31 .
To explore the molecular mechanism underlying sponge gourd fruit skin coloration, we conducted an integrated analysis of the transcriptome, sRNAome, and degradome of two sponge gourd materials with distinct fruit skin colors (i.e., green and white). The key miRNAs, genes, and the network associated with chlorophyll and flavonoid metabolism, chloroplast development, and chloroplast protection were identified. These results provide researchers with valuable information regarding fruit skin coloration and its complex effect on sponge gourd fruit quality.

Materials and methods
Plant materials. Two sponge gourd inbred lines, WS and GS, with distinct fruit skin colors, were used in this study. Mature WS fruits had a white peel with a green base, whereas mature GS fruits had a green peel and base. These materials were obtained and preserved in our team, the Specialty Vegetable Breeding Laboratory, Institute of Vegetables, Zhejiang Academy of Agricultural Sciences. Plant materials were planted at the Haining Innovation Base of Zhejiang Academy of Agricultural Sciences, under standard field management. The fruit skin samples were collected from the middle part of the fruits at the mature fruit stage, then were immediately frozen in liquid nitrogen and stored at − 80 °C. Three biological replicates were analyzed for both materials. The collection of plant material complied with the relevant institutional, national, and international guidelines and legislation. 1.0 g fruit skin samples per replicate were weighed and ground into powder in liquid nitrogen. The ground material was added to 10 mL ethanol and incubated at room temperature until the material turned completely white. The samples were passed through filter paper. The filtrate was topped up with ethanol for a final volume of 30 mL. The absorbance of the samples was measured at 649 nm (A649) and 665 nm (A665) using a spectrophotometer, with ethanol as the blank control. Three replicates of chlorophyll contents were analyzed for both GS and WS lines. The chlorophyll contents were calculated and the data were analyzed by T-test using SAS 8.0 software with p < 0.01.
Transmission electron microscopy analysis. To observe the chloroplast ultrastructure of the fruit skin, WS and GS fruit skin tissues were excised using a sterile razor blade and then immediately fixed in 2.5% (v/v) glutaraldehyde solution. The samples were subsequently rinsed with 0.1 M phosphate buffered solution (PBS) buffer, fixed in 1.0% (v/v) osmium tetroxide, rinsed with 0.1 M PBS buffer, dehydrated using graded ethanol, embedded using the Spurr kit, and polymerized at 70 °C. Ultrathin sections cut using a Leica UC6 ultramicrotome (Leica, Germany) were stained with uranyl acetate for 15 min and then with lead citrate for 5 min at room temperature. Chloroplasts were examined using the H7650 microscope (Hitachi, Tokyo, Japan).
mRNA library construction and sequencing. Total RNA was isolated and purified using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's recommended procedure. Poly-(A) RNA was purified from 1 μg total RNA using Dynabeads Oligo-(dT)25-61005 (Thermo Fisher, CA, USA), with two rounds of purification. The poly-(A) RNA was fragmented at 94 °C using the Magnesium RNA Fragmentation Module (NEB, USA). The RNA fragments were reverse transcribed to cDNA using SuperScript II Reverse Transcriptase (Invitrogen). The first-strand cDNA was then used to synthesize U-labeled second-stranded cDNA. An A-base was then added to the blunt ends of each strand to prepare them for the ligation to the index adapters.
Single-or dual-index adapters were ligated to the fragments before the size selection step was performed using AMPureXP beads. After the heat-labile UDG enzyme (NEB, USA) treatment of the U-labeled cDNA, the ligated products were amplified by PCR. The average insert size for the final cDNA library was 300 ± 50 bp. Finally, the cDNA library was sequenced (2 × 150-bp paired-end sequencing) on the Illumina NovaSeq 6000 system (LC-Bio, Hangzhou, China) following the vendor's recommended protocol.
Sequence mapping and bioinformatic analysis of mRNAs. The Cutadapt software 32 was used to remove reads with adapters. After eliminating the low-quality reads and the reads with undetermined bases, the HISAT2 software 33 was used to map the remaining reads to the Luffa cylindrica genome 34 . The mapped reads of each sample were assembled using the default parameters of StringTie 35 . The transcriptomes of all samples were merged to construct a comprehensive transcriptome using the gffcompare software (https:// github. com/ gpert ea/ gffco mpare). After generating the final transcriptome, StringTie 35 and ballgown 36 were used to estimate the expression levels of all genes, which were determined in terms of fragments per kilobase of transcript per million fragments mapped (FPKM) values 37 . The differentially expressed genes (DEGs) were identified using DESeq2, with the following criteria: fold-change > 2 or < 0.5 and p < 0.05 38 . The functions of the DEGs were determined on the basis of GO enrichment 39 and KEGG enrichment 40 analyses. Heatmap of DEGs bioinformatic analysis was performed using the OmicStudio tools at https:// www. omics tudio. cn/ tool. sRNA library construction and sequencing. TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, USA) were used to prepare sRNA sequencing libraries, which were then sequenced (1 × 50-bp single-end sequencing) on the Illumina HiSeq 2500 system (LC-Bio). sRNA sequence analysis. Raw reads were analyzed using the in-house program ACGT101-miR (LC Sciences, Houston, TX, USA) to remove adapter dimers, junk, low-complexity sequences, common RNA families (e.g., rRNAs, tRNAs, snRNAs, and snoRNAs), and repeats. Unique sequences 18-25 nt long were mapped to precursors from specific species in the miRBase 22.0 database 41 to identify known miRNAs and novel 3p-and 5p-derived miRNAs. Length variations at the 3′ and 5′ ends and one mismatch within the sequence were allowed for the alignment. The unique sequences mapped to the hairpin arms of mature miRNAs from specific species were designated as known miRNAs. The unique sequences mapped to the other arm of known precursor hairpins opposite of the annotated mature miRNA-containing arm were designated as novel 5p-or 3p-derived miRNA candidates. The remaining sequences were mapped to precursor hairpins from other selected species (i.e., specific species were excluded) in the miRBase 22.0 database via a BLAST search. The mapped pre-miRNAs were used as queries in a BLAST search of the genomes from specific species to determine their genomic loca- www.nature.com/scientificreports/ tions. These sequences were defined as known miRNAs. The unmapped sequences served as queries for a BLAST search of specific genomes. The hairpin RNA structures containing these sequences were predicted according to the flanking 120-nt sequences using the RNAfold software.

Analysis of differentially expressed miRNAs (DE-miRNAs).
The DE-miRNAs identified on the basis of normalized deep-sequencing read counts were analyzed using the T-test. The significance threshold was set at 0.05. cDNA library construction for degradome sequencing. Total RNA was isolated and purified using the TRIzol reagent (Invitrogen) following the manufacturer's recommended protocol. Poly-(A) RNA was purified from the total RNA (20 µg) using oligo-(dT) magnetic beads. The 5′ adapters were ligated to the 5′ end of the 3′ mRNA cleavage products using RNA ligase. First-strand cDNA was synthesized by reverse transcription involving a 3′-adapter random primer prior to the size selection step performed using AMPureXP beads. The cDNA was amplified by PCR. The average insert size for the final cDNA library was 200-400 bp. Finally, the cDNA library was sequenced (1 × 50-bp single-end sequencing) on the Illumina HiSeq 2500 system (LC-Bio) following the vendor's recommended procedure.
Degradome sequencing data processing and target identification. The raw data were processed to obtain sequences suitable for the subsequent analysis. The sequences were aligned with the cDNA sequences in the database of sequenced species to produce the degradome density file. The mRNA sequences of the target genes paired with the sRNA sequences were predicted using the cleavage site prediction software CleaveLand4:GSTAr. The predicted target genes corresponding to miRNAs were combined with the mRNAs in the degradome density file to identify the common mRNAs, which were designated as the miRNA target genes.  Table S1.
Subcellular localization analyses of selected proteins. The CDS sequences of ten selected proteins from GS or WS were amplified and cloned into the pFGC-eGFP plasmid using One-step Fusion Cloning Mix (Toroivd, Japan) using gene-specific primer with Bam HI cleavage site (Table S2). These recombinant plasmids were transformed into Agrobacteriumt tumefaciens GV3101 and transiently expressed in tobacco leaf cells. Images were acquired at 48 h using a Leica DMLE camera (Leica, Wetzlar, Germany).
Ethical approval. This article does not contain any studies with human participants or animals performed by any of the authors.

Results
Phenotypes, chloroplast ultrastructures, and chlorophyll contents of WS and GS. Mature fruits of white skin (WS) are presenting white fruit skin color, with green color at the fruit base, and mature fruits of green skin (GS) are presenting green fruit skin color (Fig. 1a). A transmission electron microscopy analysis revealed there were fewer chloroplasts, with fewer thylakoids per chloroplast, in fruit skin cells of WS than in GS (Fig. 1a). The chlorophyll contents of fruit skin differed significantly between GS and WS. Specifically, the chlorophyll a, chlorophyll b, and total chlorophyll contents were respectively 0.266, 0.078, and 0.344 mg/g for GS, but were 0.0268, 0.0191, and 0.0459 mg/g for WS (Fig. 1b). Identification and functional annotation of DEGs. In total, 14,797 genes were detected as expressed in at least one library ( Fig. 2a and Table S3), of which 14,080 were expressed in both GS and WS fruit skins, whereas 474 and 243 genes were specifically expressed in GS and WS fruit skins, respectively (Fig. 2a). On the basis of |log 2 (fold-change)|> 1 and p < 0.05, 4331 genes that were differentially expressed between GS and WS  (Table S3). The volcano map and heatmap for these DEGs are presented in Fig. 2b,c. The DEGs included 215 transcription factor genes, of which 132 and 83 were down-regulated and up-regulated, respectively, in WS fruit skins (relative to the expression levels in GS fruit skins) (Table S4) (Fig. 3a). The significantly enriched KEGG pathways were photosynthesis-antenna proteins (ko00196), plant hormone signal transduction (ko04075), phenylpropanoid biosynthesis (ko00940), galactose metabolism (ko00052), flavonoid biosynthesis (ko00941), carbon fixation in photosynthetic organisms (ko00710), arginine biosynthesis (ko00220), glycosaminoglycan degradation (ko00531), carotenoid biosynthesis (ko00906), monoterpenoid biosynthesis (ko00902), and sesquiterpenoid and triterpenoid biosynthesis (ko00909) (Fig. 3b).
DEGs involved in chlorophyll metabolism. The expression levels of genes involved in chlorophyll metabolism were analyzed (Fig. 4). The expression levels of chlorophyll biosynthesis-related genes, including HEMA (glutamyl-tRNA reductase 1, Maker00033993), CHLM (magnesium protoporphyrin IX methyltransferase, Maker00007651), CRD1 (magnesium protoporphyrin IX monomethyl [oxidative] ester cyclase, Maker00013112), POR (protochlorophyllide reductase, Maker00016117 and Maker00000841), and CAO (chlorophyllide a oxygenase, Maker00008808), were significantly down-regulated in WS ( Fig. 4 and Table 2). Additionally, the expression levels of genes contributing to chlorophyll degradation, such as CLH1 (chlorophyllase-1-like, Maker00038018) and two SGR genes (protein STAY-GREEN, Maker00003033 and Maker00008858), were also down-regulated in WS ( Fig. 4 and Table 2). The fact that the expression of all these genes involved in chlorophyll metabolism was down-regulated in WS implied chlorophyll biosynthesis was repressed in the fruit skins of WS sponge gourd.
Twelve DEGs related to photosynthesis were detected ( Table 2). With the exception of the up-regulated expression of Maker00005035 (photosystem II stability/assembly factor), the expression levels of other genes, including two genes (Maker00013010 and Maker00004248) encoding the oxygen-evolving enhancer protein, five genes (Maker00017083, Maker00012955, Maker00012025, Maker00017196, and Maker00019563) encoding the photosystem I reaction center subunit, one gene (Maker00014866) encoding the photosystem II reaction center W protein, one gene (Maker00013798) encoding the photosystem II stability/assembly factor, and two genes (Maker00025426 and Maker00032710) encoding the psbP domain-containing protein, were down-regulated in Table 1. Summary of the RNA-Seq data for the fruit skins of GS and WS.     (Table 2). Chloroplast protection is largely associated with peroxidase gene expression levels. In this study, 25 peroxidase genes were differentially expressed between the WS and GS fruit skins, of which the expression levels of 21 and 4 genes were respectively down-regulated and up-regulated in WS fruit skins (relative to the expression levels in GS fruit skins). These results indicated that the expression of genes involved in chloroplast development and protection was mostly repressed in the WS fruit skins, thereby limiting chlorophyll biosynthesis.
DEGs involved in flavonoid biosynthesis. Twenty DEGs (18 down-regulated and 2 up-regulated) associated with flavonoid biosynthesis were detected in this study (Table 3) Table 4). The length distribution of the sRNA sequences (18-25 nt) was analyzed (Fig. 5a). More than 85% of all sRNA sequences were 21-24 nt long in the GS and WS fruit skin libraries, with 24-nt being the most common length (Fig. 5a). The sRNA sequences were used as queries for a BLAST search of the Rfam database to obtain the noncoding RNAs (e.g., rRNAs, tRNAs, snoRNAs, snRNAs, and others). For all six libraries, rRNAs were the most abundant non-coding RNAs, followed by the tRNAs, snoRNAs, snRNAs, and others (Fig. 5b).
Identification of miRNAs and analysis of their differential expression. A total of 862 miRNAs, with 826 pre-miRNA sequences, were predicted in the analyzed libraries (Table S5). Details regarding these miR-

Network of DE-miRNAs and their differentially expressed target genes involved in flavonoid biosynthesis. Several DE-miRNAs and their differentially expressed target genes involved
in flavonoid biosynthesis were identified (Table 6). For example, cme-MIR169r-p3, cme-miR396a, cme-miR159a_R-1, and PC-5p-39803_80 targeted four genes (Maker00015252, Maker00033913, Maker00039139, and Maker00034132) encoding bHLH transcription factors. Additionally, cme-MIR319b-p5 targeted two MYB Table 5. Differentially expressed miRNAs and their differentially expressed target genes associated with chlorophyll metabolism and chloroplast development. Subcellular localization analyses of selected proteins. The subcellular localization of ten selected proteins from GS and WS were analyzed by transient expression of the green fluorescent protein (GFP) fusion proteins in tobacco leaf epidermal cells. As shown in Fig. 9, six proteins from GS, including Maker00022517, Maker00022071, Maker00003033, Maker00008808, Maker00017083 and Maker00007651, were localized in the

Discussion
Researchers have focused on fruit skin color-related traits and revealed the associated genes and potential molecular mechanisms in cucumber [10][11][12][13] , tomato [14][15][16] , and potato 24 . However, sponge gourd fruit skin coloration remains relatively unexplored. This study involved an integrated investigation of the transcriptome, sRNAome, and degradome of two sponge gourd lines that differed in terms of their fruit skin colors (i.e., white and green). The generated data revealed several DE-miRNAs and DEGs involved in chlorophyll and flavonoid metabolism or chloroplast development and protection that influence sponge gourd fruit skin coloration. The chlorophyll biosynthetic genes HEMA, CHLM, CRD1, POR, and CAO were expressed at lower levels in WS fruit skins than in GS fruit skins. Chlorophyll biosynthesis is tightly regulated by HEMA; this gene encodes a glutamyl-tRNA reductase, which is the initial enzyme of the rate-limiting step involving the synthesis of 5-aminolevulinic acid (ALA) 43 . In rice, YGL18, which encodes a magnesium protoporphyrin IX methyltransferase (CHLM), is essential for light-related chlorophyll synthesis and light intensity-associated plant growth 44 . The rice PGL gene encoding chlorophyllide a oxygenase 1 (CAO1) is mainly expressed in the chlorenchyma and is activated in the light-dependent chlorophyll synthesis process; compared with wild-type plants, pgl mutant plants have a lower chlorophyll content and a disordered thylakoid ultrastructure 45 . In the current study, the WS fruit skins had lower chlorophyll contents and altered chloroplast ultrastructures compared with the GS fruit skins. The down-regulated expression of these biosynthetic genes may help to explain the inhibited chlorophyll biosynthesis in WS fruit skins.
The expression of SGR and CLH1, which are involved in chlorophyll degradation, was also down-regulated in WS fruit skins. The SGR gene was initially identified in pea as a key regulator of chlorophyll degradation that is responsible for Mendel's green cotyledon trait 46 . Mutations to STAY-GREEN-encoding homologs are responsible for the green flesh and chlorophyll retainer phenotypes of tomato and pepper 47 . Chlorophyllase (CLH) is the first enzyme in the chlorophyll catabolic pathway 48 . Earlier research identified CLH1 as a chlorophyll dephytylase www.nature.com/scientificreports/ involved in PSII repair in Arabidopsis 49 . The observed down-regulation in SGR and CLH1 expression in WS fruit skins suggests chlorophyll degradation is also impaired in WS sponge gourd fruit skins. In this study, DEGs associated with chloroplast development, including those encoding BEL1-like homeodomain protein, homeobox knotted-1-like protein, ARF, and ARR17-like protein, were identified. In ripening tomato fruits, BEL1-LIKE HOMEODOMAIN4 influences chlorophyll accumulation, chloroplast development, cell wall metabolism, and carotenoid accumulation 50 . In the WS fruit skins, the Maker00014487 (BEL1-like homeodomain protein 7) expression level was up-regulated, in contrast to the down-regulated expression of the remaining four BEL1-like genes. The KNOX genes TKN4 and TKN2 function upstream of GOLDEN2-LIKE 2 (SlGLK2) and the related gene SlAPRR2-LIKE to influence various chloroplast developmental processes in tomato fruits 17 . In this study, six KNOX genes and one ARR17-like gene were expressed at lower levels in WS fruit skins than in GS fruit skins. Previous studies demonstrated that SlARF4 negatively regulates chlorophyll accumulation in tomato fruits 19,20 , whereas SlARF10 positively regulates chlorophyll accumulation by activating SlGLK1 expression 51 . Four ARF-encoding genes, ARF3, ARF4, ARF6, and ARF8, were more highly expressed in fruit skins of WS than in GS, whereas the opposite expression pattern was observed for the ARF19-like gene. These findings suggest BEL1-like and KNOX genes may enhance chlorophyll accumulation and chloroplast development in sponge gourd, whereas ARF3, ARF4, ARF6, and ARF8 have the opposite effects.
In present study, the miR156, miR159, miR172, miR319, and miR396 family members were mostly upregulated, whereas the miR166, miR167, miR171, miR390, and miR393 family members were mostly downregulated in WS fruit skins (Table S5). Previous investigations determined miR156 is associated with anthocyanin accumulation in the pear fruit peel and in blueberry 22,23 . Additionally, miR396, miR156, miR171, and miR319 are reportedly associated with chlorophyll metabolism in plants 23,[52][53][54] . For example, transgenic creeping bentgrass overexpressing Osa-miR396c (i.e., a rice miRNA396 gene) develops abnormally and accumulates more chlorophyll than the wild-type control 52 . The overexpression of the blueberry gene VcMIR156a in tomato enhances anthocyanin biosynthesis and chlorophyll degradation in the stem by altering pigment-associated gene expression, while also altering the chloroplast ultrastructure 23 . In Arabidopsis, the miR171-SCL module is critical for mediating GA-DELLA signaling during the coordinated regulation of chlorophyll biosynthesis 53 . Transgenic tomato plants overexpressing wild tomato (Solanum habrochaites) miRNA319d (sha-miR319d) exhibit enhanced stress tolerance and have high chlorophyll contents 54 . The miRNAs identified in the current study, including miR396, miR156, miR171, and miR319, are potential regulators of sponge gourd fruit skin coloration.
Eight DE-miRNAs targeting eight DEGs with opposite expression level and associated with chlorophyll binding and chloroplast development in sponge gourd were detected (Table 5). For example, the down-regulated miRNAs aly-miR167d-3p_2ss9CT20GA and cme-miR167d targeted the up-regulated Maker00008400 (ARF8) gene. More specifically, miR167, which is an important regulator of auxin-mediated development, reportedly targets the members of a large family of transcription factors that modulate gene expression in response to auxin (i.e., ARF4, ARF6, and ARF8) 55,56 . Earlier investigations confirmed ARF4 controls fruit chloroplast development or chlorophyll biosynthesis in tomato 19,20 . In tobacco, NtARF8 regulates ANS and DFR expression in an NtTTG2-dependent manner, thereby contributing to anthocyanin production and flower coloration 57 . However, the involvement of miR167 and its target gene ARF8 in chloroplast development or chlorophyll biosynthesis remains unknown.
Six miRNAs targeting seven DEGs with opposite expression level and associated with flavonoid synthesis in sponge gourd were identified (Table 6). For example, cme-miR396a negatively targeted bHLH143, cme-miR159a_R-1 negatively targeted bHLH51, cme-MIR319b-p5 negatively targeted MYB, cme-miR319c_R + 2_1ss20TC negatively targeted CHS2, and cme-miR396b and han-miR3630-3p_L-3 negatively targeted 4CL2. The bHLH proteins form one of the largest transcription factor families in plants. The anthocyanin biosynthesis-related SGIIIf bHLH genes have been identified in horticultural crops, including those encoding bHLH3 and bHLH64 transcription factors in apple and pear 58,59 . The R2R3-MYB transcription factor genes SlAN2-like and BrPAP1a encode regulators of anthocyanin production in tomato and turnip 60,61 . Additionally, CHS and 4CL are structural genes involved in anthocyanin biosynthesis. Accordingly, these miRNAs and their targets regulate flavonoid synthesis to determine sponge gourd fruit skin coloration.

Conclusion
This study presented an integrated analysis of the transcriptome, sRNAome, and degradome between two materials of sponge gourd with distinct fruit skin colors to elucidate the genes, miRNAs and their network regulating fruit skin coloration. The crucial genes involved in chlorophyll metabolism, chloroplast development and chloroplast protection, including HEMA, CHLM, CRD1, POR, CAO, CLH, SGR, CAB, BEL1-like, KNAT and ARF, were identified. Moreover, the miR156, miR159, miR166, miR167, miR172, miR393 and their target genes involved in chlorophyll metabolism, chloroplast development and chloroplast protection were obtained. Additionally, the miR159, miR166, miR169, miR319, miR390, miR396 and their targets CHS, 4CL, bHLH and MYB involved in flavonoid biosynthesis regulatory network were identified. These results provided the molecular mechanism of fruit skin coloration at the levels of transcriptome, sRNAome and degradome, and would lay foundation for further validation of key genes and miRNAs regulating fruit skin coloration in sponge gourd.

Data availability
Raw sequencing data of transcriptome can be accessed through the GSA of National Genomics Data Center (NGDC) (https:// ngdc. cncb. ac. cn/) with the accession number CRA005331.